2 Buffer Analysis

2.1 Environment Setup

2.1.1 Load libraries

First, load the required libraries.

library(sf)
library(tidygeocoder)
library(tmap)

2.1.2 Load Data

We will use a the shapefile of spatially enabled resource data, generated from our previous tutorial. A copy of the shapefile can also be found here, but note that all four files are required (.dbf, .prj, .shp, and .shx) to consitute a shapfile. This shapefile includes methadone clinics in Chicago as our working example.

Let’s take a look at the first few rows of the dataset. Our data includes addresses but not geographic coordinates.

# Load Resource locations (shapefile)
MetClinics <- st_read("methadoneClinics.shp")

## Reading layer `methadoneClinics' from data source `/Users/maryniakolak/code/opioid-environment-toolkit/methadoneClinics.shp' using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -87.7349 ymin: 41.68698 xmax: -87.57673 ymax: 41.96475
## CRS:            4326

head(MetClinics)

## Simple feature collection with 6 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -87.67818 ymin: 41.89268 xmax: -87.63409 ymax: 41.96475
## CRS:            4326
##   X                                                         Name                 Address    City State   Zip
## 1 1                Chicago Treatment and Counseling Center, Inc. 4453 North Broadway st. Chicago    IL 60640
## 2 2                      Sundace Methadone Treatment Center, LLC 4545 North Broadway St. Chicago    IL 60640
## 3 3 Soft Landing Interventions/DBA Symetria Recovery of Lakeview    3934 N. Lincoln Ave. Chicago    IL 60613
## 4 4                                        PDSSC - Chicago, Inc.     2260 N. Elston Ave. Chicago    IL 60614
## 5 5                          Center for Addictive Problems, Inc.        609 N. Wells St. Chicago    IL 60654
## 6 6                                Family Guidance Centers, Inc.     310 W. Chicago Ave. Chicago    IL 60654
##                                    fullAdd geo_method                   geometry
## 1 4453 North Broadway st. Chicago IL 60640        osm POINT (-87.65566 41.96321)
## 2 4545 North Broadway St. Chicago IL 60640        osm POINT (-87.65694 41.96475)
## 3    3934 N. Lincoln Ave. Chicago IL 60613     census POINT (-87.67818 41.95331)
## 4     2260 N. Elston Ave. Chicago IL 60614     census POINT (-87.67407 41.92269)
## 5        609 N. Wells St. Chicago IL 60654     census POINT (-87.63409 41.89268)
## 6     310 W. Chicago Ave. Chicago IL 60654     census POINT (-87.63636 41.89657)

Next, we load a shapefile of Chicago zip code areas. You can often find shapefiles (or KML, geojson spatial data formats) on city data portals for direct download. In our case, use the data provided here, though we will walk you through downloading boundaries directly through the Census via R in a later tutorial.

# Load Zip Code locations (shapefile)
Areas <- st_read("data/chicago_zips.shp")

## Reading layer `chicago_zips' from data source `/Users/maryniakolak/code/opioid-environment-toolkit/data/chicago_zips.shp' using driver `ESRI Shapefile'
## Simple feature collection with 85 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -88.06058 ymin: 41.58152 xmax: -87.52366 ymax: 42.06504
## CRS:            4326

head(Areas)

## Simple feature collection with 6 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -88.06058 ymin: 41.73452 xmax: -87.58209 ymax: 42.04052
## CRS:            4326
##   ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10  ALAND10 AWATER10  INTPTLAT10   INTPTLON10
## 1     60501   60501        B5   G6350          S 12532295   974360 +41.7802209 -087.8232440
## 2     60007   60007        B5   G6350          S 36493383   917560 +42.0086000 -087.9973398
## 3     60651   60651        B5   G6350          S  9052862        0 +41.9020934 -087.7408565
## 4     60652   60652        B5   G6350          S 12987857        0 +41.7479319 -087.7147951
## 5     60653   60653        B5   G6350          S  6041418  1696670 +41.8199645 -087.6059654
## 6     60654   60654        B5   G6350          S  1464813   113471 +41.8918225 -087.6383036
##                         geometry
## 1 MULTIPOLYGON (((-87.86289 4...
## 2 MULTIPOLYGON (((-88.06058 4...
## 3 MULTIPOLYGON (((-87.77559 4...
## 4 MULTIPOLYGON (((-87.74205 4...
## 5 MULTIPOLYGON (((-87.62623 4...
## 6 MULTIPOLYGON (((-87.64775 4...

2.2 Simple Overlay Map

We can plot these quickly to ensure they are overlaying correctly. If they are, our coordinate systems are working correctly.

## 1st layer (gets plotted first)
tm_shape(Areas) + tm_borders(alpha = 0.4) + 
  
  ## 2nd layer (overlay)
  tm_shape(MetClinics) + tm_dots(size = 0.1, col="red") 

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

2.3 Spatial Transformation

Check the Coordinate System/Projection for your data.

st_crs(MetClinics) 

## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

Are the coordinate systems for crime points and tracts the same?

st_crs(Areas)

## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

2.3.1 Transform CRS

To calculate buffers, we will need to conver to a different CRS that preserves distance. Trying using a search engine like Google with search terms “CRS Illinois ft”, for example, to look for a code that provides what we need. After searching, we found EPSG:3435 that uses feet for a distance metric. We’ll use that!

First, set a new CRS.

CRS.new <- st_crs("EPSG:3435")

Next, transform both datasets to your new CRS.

MetClinics.3435 <- st_transform(MetClinics, CRS.new)
Areas.3435 <- st_transform(Areas, CRS.new)

Check the CRS of both datasets again. If they are identical you’re ready to move onto the next step!

2.4 Generate Buffers

# Create 2-mile buffers for each house point
MetClinic_buffers <- st_buffer(MetClinics.3435, 10560)

Inspect the structure of the object you just created. Note that this is a new data object, represented as a polygon. Each buffer around each point is a separate entity.

2.4.1 Visualize buffers

# Map in tmap
tm_shape(Areas.3435) + tm_borders() +
  tm_shape(MetClinics.3435) + tm_dots(col = "blue") +
  tm_shape(MetClinic_buffers) + tm_borders(col = "red") 

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.

# Map Housing Buffers
tm_shape(Areas) + tm_borders() +
  tm_shape(MetClinic_buffers) + tm_fill(col = "blue", alpha = .4) + tm_borders(col = "blue") +
  tm_shape(MetClinics.3435) + tm_dots(col = "red") 

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.

2.4.2 Merge buffers

union.buffers <- st_union(MetClinic_buffers)

2.4.3 Visualize union of buffers

# Map Housing Buffers
tm_shape(Areas) +  tm_borders()+
  tm_shape(union.buffers) + tm_fill(col = "blue", alpha = .2) + tm_borders(col = "blue") +
  tm_shape(MetClinics.3435) + tm_dots(col = "red", size = 0.4) 

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(gm$shape.master_crs, "proj4string"): CRS uses proj4string, which is deprecated.

## Warning in `$.crs`(crs, "proj4string"): CRS uses proj4string, which is deprecated.